元分析所用R全部代码及详细解释
目录
1. 前期准备
1.1 R的下载及安装
1.1.1 R下载及安装
1.2 Rstudio下载安装及载入程序包
1.2.1 下载安装studio
1.2.2 新建代码文件
1.2.3 设置工作路径
1.3 数据准备
1.4 读入数据
1.4.1 方式一：代码读入
1.4.2 方式二：通过Rstudio 窗口工具栏import
1.4.3 安装元分析所用程序包
2. 初步分析
2.1 聚合的效应量overall effect size
2.1.1 计算variance
2.1.2 计算总的效应量
2.1.3 未矫正的效应量
2.1.4 Baujat plot
2.1.5 Forest plot
2.1.6 Sensitivity analysis
3. 调节变量分析
3.1 连续性调节变量continuous moderators
3.1.1 数据的中心化
3.1.2 一个调节变量
3.1.3 多个连续性调节变量、
3.2 类别变量Categorical moderators
3.2.1 Dummy coding
3.3 Subgroup analysis
4. 元分析高阶篇：
4.1 多水平元分析Three level meta-analysis
4.2 跨文化元分析
4.3 结构方程模型：MetaSEM
1. 前期准备
注：这里将从软件安装到数据准备逐步说明，不假定需要有任何R使用经验。下载安装的可以略过此步骤。
1.1 R的下载及安装
1.1.1 R下载及安装
https://www.r-project.org/
下载后可以安装到本地文件夹。
1.2 Rstudio下载安装及载入程序包
1.2.1 下载安装studio
建议安装R studio：数据可视化、代码编写更容易。本篇文章将通过Rstudio进行解释说明。
任意浏览器搜索，下载、安装，打开Rstudio。
1.2.2 新建代码文件
File-new file-R Script, 如下图，或者通过快捷键Ctril+Shift+N
1.2.3 设置工作路径
在新的代码文件，将工作路径设置为数据存放的文件夹。
setwd("C:/Users/…/Data analysis")
可以先找到数据文件夹，然后复制路径。。。
F:\Phd In UtrechtUniversity
\p.li\Mymeta-analysis Leadership and engagement\Dataanalysis
注意将\换成/，然后加上双引号
setwd("F:/Phd In Utrecht University/p.li/Mymeta-analysis Leadership and engagement/Data analysis")
写好之后，选中，然后点击run就会运行
提示：为了方便最好将数据文件和新建的代码放在相同文件夹。
1.3 数据准备
R可以直接读入EXCEL数据以及SPSS 数据，建议将Excel数据存为txt格式。
1.4 读入数据
1.4.1 方式一：代码读入
##1.2 read data
dat <-read.delim(file = "llsengagement190.txt", header = TRUE)
主：本文蓝色字体均表示R代码。dat, 为了后期方便分析，将数据文件命名为dat, file 后面的即为你所命名的文件，格式为dat,sav, txt, 读入方式会有所不同。如csv
#"dataset.csv", which is in the commadelimited format
dataset <-read.csv("dataset.csv")
header = True, 是用来表示读取第一行变量名命名，否则将不会读取第一行变量名，变量名会由R自动生成。#在R里用来注释说明，#后面的语句R在运行时会自动忽略，且会显示为绿色，如下图
小tip: 如果想要将数据拆分为新的数据文件，可以通过
#usingsplit data method
addmargins(table(dat$leaadershipstyles))
dat1= split(dat, dat$leaadershipstyles)
#charismatic
dat=dat1$`11`
根据Leadership styles 将数据拆分，仅针对charismatic leadership 进行单独分析，命名为dat
1.4.2 方式二：通过Rstudio 窗口工具栏import
可以自动读入txt, exxcel, spss, stata的数据文件。
1.4.3 安装元分析所用程序包
元分析有很多程序包，建议安装metafor和 metaSEM
install.packages(c("metaSEM","metafor"))
然后安装后需要读取，以后使用的时候都需要调用程序包用library命令实现
library("metaSEM")
library(metafor)
2. 初步分析
2.1 聚合的效应量overall effect size
2.1.1 计算variance
注：本文演示的效应量是相关系数r，其它的实验或干预研究略有不同。
首先要将r转换为Z值，然后根据样本量计算variance,
#1.3 ## # transform r to Z and calculate thecorresponding sample variances.
dat <-escalc(measure="ZCOR", ri=ri, ni=ni, data=dat, slab=paste(authors,year, sep=", "))
然后查看数据就会发现多了两列，yi, vi
View(dat)
就可以查看数据，见下图。
2.1.2 计算总的效应量
#Fit a random model
####3.2Estimating summary effect##
res <-rma(yi, vi, data=dat)
res
predict(res,digits=3, transf=transf.ztor)
confint(res)
随后通过predict命令将Z值转换为r，得到数据结果
对于数据结果解读，请参考。。。r practical tutorial on conducting Meta-analysisin R
(Del Re, 2015)以及(Quintana, 2015)。
2.1.3 未矫正的效应量
#using uncorrected r
overal2 <-rma(ri, vi, data=dat)
overal2
2.1.4 Baujat plot
Q-statistics以及I2可以告诉我们Heterogeneity是否显著，但并不能告诉是哪些因素造成了影响， they donot provide information on which studies may be influencing to overallheterogeneity. If there is evidence of overall heterogeneity, construction of aBajaut plot can illustrate studies that are contribute to overall heterogeneityand the overall result. Study IDs are used to identify studies
b_res <- rma(yi, vi, data=dat,slab=studyid) # # Newmeta-analysis with study ID identifier
# The next command will plot a Baujat plot.
baujat(b_res)
越是右上角表示影响越大，可以看出30,25,23等对于heterogeneity 造成了显著的影响。查看这几个研究发现相关系数，r = .679, .806, -.57 偏高及出现了负值。
但是是不是确实显著，并不能通过直观来获得，
# Studies that fall to the top right quadrant of the Baujat plot contribute most to both thesefactors. A closer look the characteristics of these studies may revealmoderating variables that may contribute to heterogeneity
# A set ofdiagnostics are also available to identify potential outliers and influentialcases.
inf <-influence(res)
print(inf)
plot(inf)
红色的表示确实有显著影响，可以看出30为outlier，可以做|sensitivity analysis删掉outliers查看结果是否有显著差异，另外也有人根据forest plot查看研究的效应量是否在95%CI
2.1.5 Forest plot
## cumulativemeta-analysis (in the order of publication year)
forest(res,xlim=c(-4,2), atransf=transf.ztor,
at=transf.rtoz(c(-.4,-.2,0,.2,.4,.6)),digits=c(2,3), cex=.75, font =2)
### add column headings to theplot
text(-3.68, 36, "Authors andYear", font=2, cex= .8)
#36 is the height, if you haveless study like 15 so basically you can use 17
text(-3.68,"Authors and Year", pos=4, font=2, cex= .8)
text ( 2, 36,"[95% CI]", pos=2)
注意：标黄数值需要根据自己的样本量进行调节，否则会出现看不到标题的情况（见下图）。
2.1.6 Sensitivity analysis
根据forest plot以及前面的分析，删除效应量不在95%区间研究。
在R中可以删除一行或者多行来实现
#authentic
dat_bias <-dat[-c(23, 25, 30), ]
创建了新的数据文件dat_bias 删掉了三个数据23,25,30,然后再次运行计算总的效应量。
#unbiasedeffect size senssensitivity check the final step after moderation analysiss
res.b <- rma(yi, vi, data=dat_bias)
res.b
predict(res.b, digits=3,transf=transf.ztor)
confint(res.b)
2.1.7 Funnel plot
### funnel plot 用来检验publication bias 是否对称
funnel(res, xlab = "Correlation coefficient")
#Tests for bias
regtest(res)
ranktest(res)
3. 调节变量分析
3.1 连续性调节变量continuous moderators
3.1.1 数据的中心化
pds.c<-scale(dat $pds, center = TRUE, scale = FALSE)[,]
以power distance 为例，pds.c为创建的新的变量名，dat $pds 选择数据中原来的pds变量，将pds中心话，也就是每个数据减去其平均值。
3.1.2 一个调节变量
#4.3 power distance
res.modpds<-rma(yi, vi, mods = ~ pds.c, data=dat_bias)
res.modpds
3.1.3 多个连续性调节变量、
只需要modes = ~ 加上即可
res.modculture<- rma(yi, vi, mods = ~ pds.c + pori.c +fuori., data=dat_bias)
res.modculture
以上为三个Moderators
3.1.3 交互作用
使用乘号*即可
#4.3 power distance
res.modpds2<- rma(yi, vi, mods = ~ pds.c*ctlsenes.c, data=dat_bias)
res.modpds2
3.2 类别变量Categorical moderators
3.2.1 Dummy coding
例子：以研究设计：横断数据/纵向数据/日记研究三种类别
##6.3 study design
首先查看每个类别的数量，如下
addmargins(table(dat_bias$design))
1 2 Sum
28 2 30
发现只有两种，创建两个类别变量，
crossec<- ifelse(dat $design=="1", yes=1, no=0)
longit<- ifelse(dat $design=="2", yes=1, no=0)
summary( Model63 <- meta(y=yi, v=vi,intercept.constraint=0,
+ x=cbind(crossec,longit),
+ data=dat_bias,
+ model.name="Model 64") )
注意这里用到了metaSEM程序包，通过限定intercept 为0，即可以得到两种类别下，效应量。如果出错，或许是yi vi数据类型不对，尝试用以下代码，
dat$yi<-as.numeric(dat$yi)
dat$vi<-as.numeric(dat$vi)
有时候会数据类型成了factor,metaSEM就会出错，将数据转换为numeric 数据，或者用metafor ：
res.modpds<- rma(yi, vi, mods = ~ pds.c -, data=dat_bias)
res.modpds
添加减号，就可以删除intercept
3.3 Subgroup analysis
类似与categorical moderator.
4. 元分析高阶篇：
4.1 多水平元分析Three level meta-analysis
参考文献FittingThree-Level Meta-Analytic Models in R: A Step-by-Step Tutorial Mark (Assink & Wibbelink, 2016)。这篇文章讲的非常详细。我就不再粘贴代码了。简单说明何种情况下做\three level, 就是数据出现嵌套，一种情况在发展心理学很多研究都是一个研究通常包含了多个study，这个时候就是一种嵌套数据。或者如有人用到国家变量如GDP，国家下面有多个研究也是一种嵌套。
Assink, M., & Wibbelink, C.J. M. (2016). Fitting three-level meta-analytic models in R: A step-by-steptutorial. The Quantitative Methods for Psychology, 12(3),154–174.
4.2 跨文化元分析
很多关于工业组织及管理学领域的/领导之类的研究都可以做跨文化的元分析。常见的有两种方式：将将国家编码为categorical variable,然后比较不同区域的差异；另外一种是编码每个国家在不同文化维度上面的分数，如hofstede 研究。在网站上都可以查到很多关于每个国家在各个维度上的分数，然后拿这些分数作为连续性调节变量。就可以检验文化的影响。其实也就是实用二手数据。
https://www.hofstede-insights.com/country/china/
还有一个是关于 Culture, Leadership, and Organizations: The Globe Study of 62 Societies
类似的研究有已经发表在顶级期刊，JAP，如下
Liu, D., Jiang, K., Shalley, C.E., Keem, S., & Zhou, J. (2016). Motivational mechanisms of employeecreativity: A meta-analytic examination and theoretical extension of thecreativity literature. Organizational Behavior and Human Decision Processes,137, 236–263.
Rabl, T.,Jayasinghe, M., Gerhart, B., & Kühlmann, T. M. (2014). A meta-analysis of countrydifferences in the high-performance work system-business performancerelationship: The roles of national culture and managerial discretion. Journalof Applied Psychology, 99(6),
Rockstuhl, T., Dulebohn, J. H.,Ang, S., & Shore, L. M. (2012). Leader-member exchange (LMX) and culture: Ameta-analysis of correlates of LMX across 23 countries. Journal of AppliedPsychology, 97(6), 1097–1130.
Taras, V., Kirkman, B. L., & Steel, P. (2010). Examiningthe Impact of Culture’s Consequences: A Three-Decade, Multilevel, Meta-AnalyticReview of Hofstede’s Cultural Value Dimensions. Journal of AppliedPsychology, 95(3), 405–439.
4.3 结构方程模型：MetaSEM
这个说实在相当比较复杂。我自己还在摸索中。推荐的文献两个
AND META ANALYSIS Meta-Analytic StructuralEquation Modelling (Jak, n.d.)。
Meta-Analysis: A Structural Equation ModelingApproach (Cheung, 2015)
如果时间有限可以读下第一篇，92页，相对友好，从数据准备到后期分析都有逐步解释。Suzanne Jak 是Cheung的学生，当然后者比较全面，不仅中介调节等都有讲，也有实例R代码。
最后可以通过多种方式获得这篇文章用到的代码及数据：
1：自己输入或者copy
2：转发至朋友圈,
3：点击文章下方like the author.发送邮件
：直接发送邮件。
元分析从去年搞到今年各种幺蛾子，然而仍然还在进行中。。。。R确实非常实用，若有错误欢迎指正。